A UNIFORM INF-SUP CONDITION WITH APPLICATIONS TO 

PRECONDITIONING 

KENT-ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER 



Abstract. A uniform inf-sup condition related to a parameter dependent Stokes problem 
is established. Such conditions are intimately connected to the construction of uniform 
preconditioners for the problem, i.e., preconditioners which behave uniformly well with 
respect to variations in the model parameter as well as the discretization parameter. For 
the present model, similar results have been derived before, but only by utilizing extra 
regularity ensured by convexity of the domain. The purpose of this paper is to remove 
this artificial assumption. As a byproduct of our analysis, in the two dimensional case we 
also construct a new projection operator for the Taylor-Hood element which is uniformly 
bounded in L 2 and commutes with the divergence operator. This construction is based on 
a tight connection between a subspace of the Taylor-Hood velocity space and the lowest 
order Nedelec edge element. 



I. Introduction 

The purpose of this paper is to discuss preconditioners for finite element discretizations 
of a singular perturbation problem related to the linear Stokes problem. More precisely, let 
Q C W 1 be a bounded Lipschitz domain and e G (0, 1] a real parameter. We will consider 
singular perturbation problems of the form 

(I — e 2 A)u — gradp = / in Q, 
(1-1) div-u — g in f2, 

u — on <9f2, 

where the unknowns u and p are a vector field and a scalar field, respectively. For each fixed 
positive value of the perturbation parameter e the problem behaves like the Stokes system, 
but formally the system approaches a so-called mixed formulation of a scalar Laplace equa- 
tion as this parameter tends to zero. In physical terms this means that we are studying fluid 
flow in regimes ranging from linear Stokes flow to porous medium flow. Another motivation 
for studying preconditioners of these systems is that they frequently arises as subsystems in 
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time stepping schemes for time dependent Stokes and Navier-Stokes systems, cf. for example 

El ESI HBlHg. 

The phrase uniform preconditioners for parameter dependent problems like the ones we 
discuss here, refers to the ambition to construct preconditioners such that the preconditioned 
systems have condition numbers which are bounded uniformly with respect to the perturba- 
tion parameter e and the discretization. Such results have been obtained for the system( Jl~T] ) 
in several of the studies mentioned above, but a necessary assumption in all the studies so 
far has been a convexity assumption on the domain, cf. [IT]. However, below in Section [3] we 
will present a numerical example which clearly indicates that this assumption should not be 
necessary. Thereafter, we will give a theoretical justification for this claim. The basic tool 
for achieving this is to introduce the Bogovskii operator, cf. [Bj, as the proper right inverse 
of the divergence operator in the continuous case. 



The construction of uniform preconditioners for discretizations of systems of the form ( 1.1 ), 
is intimately connection to the well-posedness properties of the continuous system, and the 
stability of the discretization. In fact, if we obtain appropriate e-independent bounds on 
the solution operator, then the basic structure of a uniform preconditioner for the contin- 
uous system is an immediate consequence. Furthermore, under the assumption of proper 
stability properties of the discretizations, the basic structure of uniform preconditioners for 
the discrete system also follows. We refer to [15] and references given there for a discussion 
of these issues. The main tool for analyzing the well-posedness properties of saddle-point 



problems of the form (1.1) is the Brezzi conditions, cf. [HE]. In particular, the desired 



uniform bounds on the solution operator is closely tied to a uniform inf-sup condition of the 



form (2.4) stated below. Furthermore, the verification of such uniform conditions are closely 
tied to the construction of uniformly bounded projection operators which properly commute 
with the divergence operator. In the present case, these projection operators have to be 
bounded both in L 2 and in H 1 . In Section [4] we will construct such operators in the case of 
the Mini element and the Taylor-Hood element, where the latter construction is restricted 
to quasi-uniform meshes in two space dimensions. 



2. Preliminaries 



To state the proper uniform inf-sup condition for the system (1.1) we will need some 
notation. If X is a Hilbert space, then || ■ \\x denotes its norm. We will use H m = H m (Q) 
to denote the Sobolev space of functions on Q with m derivatives in L 2 = L 2 (Q). The 
corresponding spaces for vector fields are denoted H m (Q;W l ) and L 2 (f2;M n ). Furthermore, 
(•, •) is used to denote the inner-products in both L 2 (fl) and L 2 (Q;M. n ), and it will also 
denote various duality pairings obtained by extending these inner-products. In general, we 
will use H™ to denote the closure in H m of the space of smooth functions with compact 
support in Q, and the dual space of H™ with respect to the L 2 inner product by H~ m . 
Furthermore, L\ will denote the space of L 2 functions with mean value zero. We will use 
£(X, Y) to denote the space of bounded linear operators mapping elements of X to Y, and 
if Y = X we simply write C(X) instead of C(X, X). 
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If X and Y are Hilbert spaces, both continuously contained in some larger Hilbert space, 
then the intersection X R Y and the sum X + Y are both Hilbert spaces with norms given 
by 

IMIxny = IMIx + IMIr and||z||^ +y = inf (\\x\\ 2 x + ||y||y). 

z=x+y 

Furthermore, if XC\Y are dense in both the Hilbert spaces X and Y then (XdY)* = X* + Y* 
and {X + Y)* = X* n F*, cf. 0. 





The system (1.1) admits the following weak formulation: 
Find (u,p) G H^(Q;R n ) x L 2 {Q) such that 

(w,^)+e 2 (^, J Dt;) + (p,divt;> = (f,v), v e H l (Q; R n ), 

(div u,q) =(g,q), ? e ^o( fi ) 

for given data / and g. Here Dv denotes the gradient of the vector field v. More compactly, 
we can write this system in the form 

(2.2) A e = rj , where A e = 

For each fixed positive e the coefficient operator A t is an isomorphism mapping X = 
#1(0; R n )xL 2 Q (Q) ontoX* = ff-^Q; M")xL2(fi). However, the operator norm 
will blow up as e tends to zero. 

To obtain a proper uniform bound on the operator norm for the solution operator A~ 1 we 
are forced to introduce e dependent spaces and norms. We define the spaces X e and X* by 

X t = (L 2 (leH*)(tt;R n ) x ((H 1 n L 2 ) + e' 1 L 2 )(Q) 

and 

X e * = (L 2 + e' 1 H'^in-W 1 ) x (fl^neLgXfi). 
Here ifj' 1 D Lq corresponds to the dual space of i/ 1 fl Lq. Note that the space X e is equal 
to X as a set, but the norm approaches the L 2 -norm as e tends to zero. 

Our strategy is to use the Brezzi conditions |H E] to claim that the operator norms 

(2.3) \\A t \\ C (Xt,x*) and \\A; l \\ C {x:,x t ) are bounded independently of e . 

In fact, the only nontrivial condition for obtaining this is that we need to verify the uniform 
inf-sup condition 

(2.4) sup ( d j v ^g) > a \\ q \\ m+e _ lL2) qeL 2 (Q), 

where the positive constant a is independent of e G (0, 1]. Of course, if e > is fixed, and a 
is allowed to depend on e, then this is just equivalent to the standard inf-sup condition for 
the stationary Stokes problem. 



As explained, for example in [T5], the mapping property (2.3) implies that the "Riesz 
operator" B e , mapping X* isometrically to X e , is a uniform preconditioner for the operator 
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A e . More precisely, up to equivalence of norms the operator B e can be identified as the block 
diagonal and positive definite operator B t : X* — > X e , given by 



(2.5) B t 




A)- 1 + e 2 / 



This means that the preconditioned coefficient operator B e A t is a uniformly bounded family 
of operators on the spaces X e , with uniformly bounded inverses. Therefore, the precondi- 
tioned system 

B e A t 

can, in theory, be solved by a standard iterative method like a Krylov space method, with 
a uniformly bounded convergence rate. We refer to [12] for more details. Of course, for 
practical computations we are really interested in the corresponding discrete problems. This 
will be further discussed in Section S] below. 




3. The uniform inf-sup condition 



The rest of this paper is devoted to verification of the uniform inf-sup condition (2.4), and 
its proper discrete analogs. We start this discussion by considering the standard stationary 
Stokes problem given by: 

Find (u,p) G Hq(Q;R u ) x Lg(fi) such that 

(Du,Dv) + (p,dwv) =(f,v), veH^(Q;R n ), 
(divw,g) =(g,q), 9 6 ^o( fi )- 

where (/, g) G //" _1 (fi; R n ) x Lq(Q). The unique solution of this problem satisfies the estimate 
(3-2) Nlffi + IMU» <c(\\f\\ H -i + \\g\\^), 

cf. [13j. Furthermore, if the domain Q is convex, / G L 2 (Q]R n ), and g = then u e 
H 2 n Hq(Q; W 1 ), p G H 1 fl Lq(Q), and an improved estimate of the form 

(3-3) ||«||fla + \\p\\m < c H/IIra, 

holds ([9]). 



We define R G £(H R n ), Ll(Q)) to be the solution of operator of the system (3.1), 
with g = 0, given by / i— > p = Rf. Hence, if the domain Q is convex, this operator will also 
be a bounded map of L 2 (fi; R n ) into H 1 n L 2 (Q). Furthermore, let S G C{L 2 (Q), H^(Q; R n ) 
denote the corresponding solution operator, defined by (3.1 ) with / = 0, given by g h-> u = 
Sg. Then S is a right inverse of the divergence operator, and the operator S is the adjoint 
of R, since 



(f,Sg) = (Rf,dwSg) + (Du,DSg) = (Rf,g) + (p,divu) = (Rf,g). 
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Here u and p are components of the solutions of (3.1 ) with data (/, 0) and (0, g), respectively. 
As a consequence of the improved estimate (3.3 ), we can therefore conclude that if the domain 
Q is convex then S can be extended to an operator in L 2 (Q; IR™)). In other words, 

in the convex case we have 

(3.4) S G C(Lq, Hq) n C{Hq X , L 2 ), and divSg = g. 

However, the existence of such a right inverse of the divergence operator implies that the 
uniform inf-sup condition holds, since for any q G Lq(Q) we have 

(g,q) (div Sg,q) {divv,q) 



\q\\w+e-iv>= sup " " <c sup 



< c sup 



On the other hand, if the domain Vt is not convex, then the estimate (3.3) is not valid, and as 
a consequence, the operator S cannot be extended to an operator in £(if^" 1 (f2), L 2 (Q; W 1 )). 
Therefore, the proof of the uniform inf-sup condition (2.4) outlined above breaks down in 
the nonconvex case. 



3.1. General Lipschitz domains. The main purpose of this paper is to show that the 
problems encountered above for nonconvex domains are just technical problems which can 
be overcome. As a consequence, preconditioners of the form B e given by (2.5) will still 
behave as a uniform preconditioner in the nonconvex case. To convince the reader that this 
is indeed a reasonable hypothesis we will first present a numerical experiment. We consider 
the problem (1.1) on three two dimensional domains, referred to as Q2 and f^. Here fl\ 
is the unit square, f2 2 is the L-shaped domain obtained by cutting out the an upper right 
subsquare from fli, while ^3 is the slit domain where a slit of length a half is removed from 
from Qi, cf. Figure [TJ Hence, only Qi is a convex domain. The corresponding problems 



Figure 1. The domains Hi, f^, and f2 3 



(1.1 ) were discretized by the standard Taylor-Hood element on a uniform triangular grid to 



obtain a discrete anolog of this system (2.2) on the form 



A 



Uh 




Here the parameter h indicates the mesh size. We have computed the condition numbers 
of the operator B e ^A ei u f° r different values of e and h for the three domains. The operator 



B e> h is given as the corresponding discrete version of (2.5), i.e., exact inverses of the discrete 
elliptic operators appearing in (2.5) are used. Hence, in the notation of [T3] a canonical 



preconditioner is applied. The results are given in Table [T] below. 
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Table 1. Condition numbers for the operators B^hA^h 



These results indicate clearly that the condition numbers of the operators B^A^^ are not 
dramatically effected by lack of convexity of the domains. Actually, we will show below that 
these condition numbers are indeed uniformly bounded both with respect to the perturbation 
parameter e and the discretization parameter h. 



We will now return to a verification of the inf-sup condition (2.4) for general Lipschitz 



domains. The problem we encountered above in the nonconvex case is caused by the lack of 
regularity of the solution operator for Stokes problem on general Lipschitz domains. However, 



to establish (2.4) we are not restricted to such solution operators. It should be clear from 



the discussion above that if we can find any operator S satisfying condition (3.4), then (2.4) 
will hold. A proper operator which satisfy these conditions is the Bogovskii operator, see for 
example [HI section III. 3], or [8, 12]. On a domain Q, which is star shaped with respect to 



an open ball B, this operator is explicitly given as an integral operator on the form 



Sg(x) 



g{y)K(x - y, y) dy, where K(z, y) 



6(y + r—)r 

\z\ \ z 



n-1 



dr. 



Here 6 G C^°( 



with 



supp^ C B, and / 9(x) dx = 1 



This operator is a right inverse of the divergence operator, and it has exactly the desired 



mapping properties given by (3.4), cf. [HI [12]. Furthermore, the definition of the right inverse 
S can also be extended to general bounded Lipschitz domains, by using the fact that such 
domains can be written as a finite union of star shaped domains. The constructed operator 



will again satisfy the properties given by (3.4). We refer to [HJ section III. 3], [121 section 2], 
and [HI section 4.3] for more details. We can therefore conclude our discussion so far with 
the following theorem. 
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Theorem 3.1. Assume that is a bounded Lipshitz domain. Then the uniform inf-sup 



condition (2.4) holds. 

4. Preconditioning the discrete coefficient operator 



The purpose of this final section is to show discrete variants of Theorem 3.1 for various 



finite element discretizations of the problem (1.1). More precisely, we will consider finite 



element discretizations of the system (1.1) of the form: 

Find (uh,Ph) 6 K x Q h such that 
, 41 v (uh,v) + e 2 (Du h ,Dv) + (p h , divv) = (f,v), v E V h , 

(div Uh,q) =(g,q), qEQ h . 

Here Vh and Qh are finite element spaces such that Vh x Qh C Hq (O; W 1 ) x Lg(O), and h is 
the discretization parameter. Alternatively, these problems can be written on the form 



A 



C.ll 



Uh 
.Ph. 




where the coefficient operator A e ,h is acting on elements of Vh X Qh- In the examples below we 
will, for simplicity, only consider discretizations where the finite element space Vh x Qh C X e 
for all e in the closed interval [0,1]. This implies that also the pressure space Qh is a subspace 
of H 1 . The proper discrete uniform inf-sup conditions we shall establish will be of the form 

(divv , q) 

(4.2) sup —r > a\\q\\ H i+e-iL2,h, QEQh, 

where the positive constant a is independent of both e and h. Here the discrete norm 
|| ■ Hfl-i+e-i L 2 ,h is defined as 

IMlWe-i inf n (lki||^+e- 2 ||g 2 ||i 2 ), qEQh- 

9=91+92 



The technique we will use to establish the discrete inf-sup condition (4.2) is in principle 



rather standard. We will just rely on the corresponding continuous condition (2.4) and 



a bounded projection operator into the velocity space Vh- The key property is that the 



projection operator Flh commutes properly with the divergence operator, cf. (4.4) below, 
and that it is uniformly bounded in the proper operator norm. Such projection operators 
are frequently referred to as Fortin operators. 

We will restrict the discussion below to two key examples, the Mini element and the 
Taylor-Hood element. For both these examples we will construct interpolation operators 
Uh '■ L 2 (Q;M. n ) Vh which are uniformly bounded, with respect to h, in both L 2 and Hq. 
Therefore, these operators will be uniformly bounded operators in C((L 2 fl 6Hq)(Q; R n )), 
i.e., we have 

(4.3) \\nh\\c((L 2 neH^) is bounded independently of e and h. 
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It is easy to see that this is equivalent to the requirement that Uh is uniformly bounded with 
respect to h in C(L 2 ) and C(Hq). Furthermore, the operators Ilh will satisfy a commuting 
relation of the form 

(4.4) (div n h v,q) = (divv,q) v G H*(Q;R n ), q G Q h . 



As a consequence, the discrete uniform inf-sup condition (4.2) will follow from the corre- 



sponding condition (2.4) in the continuous case. We recall from Theorem 3.1 that (2.4) 



holds for any bounded Lipschitz domain, and without any convexity assumption, and as a 



consequence of the analysis below the discrete condition (4.2) will also hold without any 
convexity assumptions. The key ingredient in the analysis below is the construction of a 
uniformly bounded interpolation operator Ilh- In the case of the Mini element the con- 
struction we will present is rather standard, and resembles the presentation already done 
in pQ, where this element was originally proposed (cf. also [5j Chapter VI]). However, for 
the Taylor-Hood element the direct construction of a bounded, commuting interpolation 
operator is not obvious. In fact, most of the stability proofs found in the literature for this 
discretization typically uses an alternative approach, cf. for example [5J Section VI. 6] and 
the discussion given in the introduction of [101 . An exception is pU], where a projection sat- 



isfying (4.4) is constructed. However, this operator is not bounded in L . Below we propose 



a new construction of projection operators satisfying (4.4) by utilizing a technique for the 



Taylor-Hood method which is similar to the construction for Mini element presented below. 



The new projection operator will be bounded in both I 2 and H 1 , and hence it satisfies (4.3) 
This analysis is restricted to quasi-uniform meshes in two space dimensions. 



In the present case, the discrete inf-sup condition (4.2) will imply uniform stability of the 
discretization in the proper norms introduced above. As a consequence, we are able to derive 
preconditioners B e> h, such that the condition numbers of the corresponding operators B e ^-A-e,h 
are bounded uniformly with respect to the perturbation parameter e and the discretization 



parameter h. The operator B e> h can be taken as a block diagonal operator of the form (2.5), 
but where the elliptic operators are replaced by the corresponding discrete analogs. In fact, 
to obtain an efficient preconditioner the inverses of the elliptic operators which appear should 
be replaced by corresponding elliptic preconditioners, constructed for example by a standard 
multigrid procedure. We refer to [15], see in particular Section 5 of that paper, for a discussion 
on the relation between stability estimates and the construction of uniform preconditioners. 
In particular, the results for the Taylor-Hood method presented below explains the uniform 
behavior of the preconditioner B e ^ observed in the numerical experiment reported in Table [T] 
above. 



4.1. The discrete inf-sup condition. The rest of the paper is devoted to the construction 
of proper interpolation operators Ilh for the Mini element and the Taylor-Hood element, i.e, 
we will construct interpolation operators II h : L 2 (Q;M. n ) — > Vh such that (4.3) and (4.4) 
holds. We will assume that the domain Q is a polyhedral domain which is triangulated by 
a family of shape regular, simplicial meshes {Th} indexed by decreasing values of the mesh 
parameter h = max.TeTh ^t- Here is the diameter of the simplex T. We recall that the 
mesh is shape regular if the there exist a positive constant 70 such that for all values of the 
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mesh parameter h 

/4<7o|T|, TeT h . 

Here \T\ denotes the volume of T. 

4.1.1. The Mini element. We recall that for this element the velocity space, Vh, consists 
of linear combinations of continuous piecewise linear vector fields and local bubbles. More 
precisely, v G Vh if and only if 

V = v 1 + Ct&T, 
TeT h 

where v 1 is a continuous piecewise linear vector field, crp G W 1 , and b T G V n +i(T) is the 
bubble function with respect to T, i.e. the unique polynomial of degree n + 1 which vanish 
on dT and with f T bTdx = 1. The pressure space Qh is the standard space of continuous 
piecewise linear scalar fields. 

In order to define the operator Tlh we will utilize the fact that the space Vh can be 
decomposed into two subspaces, V h b , consisting of all functions which are identical to zero 
on all element boundaries, i.e. V^ is the span of the bubble functions, and V h x consisting of 
continuous piecewise linear vector fields. Let Tl\ : L 2 (f2; IR n ) — ► V h b be defined by, 

(n b v,z) = (v,z) Vzez h , 

where Zh denotes the space of piecewise constants vector fields. Clearly this uniquely deter- 
mines Fl\. Furthermore, a scaling argument, utilizing equivalence of norms, shows that the 
local operators JJ^ are uniformly bounded, with respect to h, in L 2 (Q;M. n ). 



The operator H\ will satisfy property (4.4) since for all v G HQ(Q;K. n ) and q G Qh, we 
have 

(4.5) (div n b h v,q) = -{n b h v, gradg) = -(w,gradg) = (divv,q), 

where we have used that gradQ^ C Z h - 

The desired operator Flh will be of the form 

n h = n b (i - R h ) + R h , 

where R h : L 2 (fl; R n ) — >■ V h x will be specified below. Note that 

i-n h = (i-n b h )(i-R h ), 

and therefore 

(div(J - n h )v,q) = (div(J - n b )(I - R h )v,q)) = 



for all q G Qh- Hence, the operator n h satisfies (4.4). 



We will take Rh to be the Clement interpolant onto piecewise linear vector fields, cf. [7]. 
Hence, in particular, the operator Rh is local, it preserves constants, and it is stable in L 2 
and Hq. More precisely, we have for any T G Th that 

(4.6) ||(I - R h )v\\ HHT) < ch k n J\\v\\ Hk{nT) , < j < k < 1, 

where the constant c is independent of h and v. Here Qt denote the macroelement consisting 
of T and all elements T' G Th such that T fl V ^ 0, and hn T = ma.XTer h ,Tcn h h?- It 
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Table 2. Condition numbers for the operators 
Mini element. 



E f hA f h discretized with the 



also follows from the shape regularity of the family {Th} that the covering {Qt}t£T h has a 
bounded overlap. Therefore, it follows from (4.6) and the L 2 boundedness of Tl\ that Tlh 
is uniformly bounded in £(L 2 (Q; W 1 )). Furthermore, by combining (4.6) with a standard 
inverse estimate for polynomials we have for any T e Th that 

||/7hu||jfi(T) < \\n b h (I - R h )v\\ H i (T) + \\R h v\\ H i {T) 



<c{h^ l \\n b h {I-R h )v\\ L 2 (T) + 



<c{h- T l \\{I 
< c{h^ l h nT 



Pllifi(T), 



Rh)V\\tf(T) + 

l)lkllifi(n T ) 

< c\\v\\m(n T ), 

, is uniformly bounded by shape regularity. This implies that 
il/j is uniformly bounded in C(Hq (O; IR n )). We have therefore verified (4.3). Together with 
(4.4) this implies (4.2). In Table [2] below we present results for the Mini element which are 



where we have used that h^h 



completely parallel to results for the Taylor-Hood element presented in Table [T] above. As we 
can see, by comparing the results of the two tables, the effect of the different discretizations 
seems to minor, as long as the mesh is the same. 



4.1.2. The Taylor-Hood element. Next we will consider the classical Taylor-Hood element. 
We will restrict the discussion to two space dimensions, and we will assume that the family of 
meshes {Th} is quasi-uniform. More precisely, we assume that there is a mesh independent 
constant 71 > such that 

(4.7) h T >ixK Te T h , 

where we recall that h = maxy hr- For the Taylor-Hood element the velocity space, Vh, 
consists of continuous piecewise quadratic vector fields, and as for the Mini element above 
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Qh is the standard space of continuous piecewise linear scalar fields. Note that if we have 
established the discrete inf-sup condition for a pair of spaces [V h ~,Qh), where V h ~ is a 
subspace of Vh, then this condition will also hold for the pair (Vh, Qh)- This observation will 
be utilized here. 

For technical reasons we will assume in the rest of this section that any T G Th has 
at most one edge in dfl. Such an assumption is frequently made for convenience when the 
Taylor-Hood element is analyzed, cf. for example [U Proposition 6.1], since most approaches 
requires a special construction near the boundary. On the other hand, this assumption will 



not hold for many simple triangulations. Therefore, in Section 4.1.3 below we will refine our 
analysis, and, as a consequence, this assumption will be relaxed. 

We let V h ~ C Vh C Hq(Q;R 2 ) be the space of piecewise quadratic vector fields which has 
the property that on each edge of the mesh the normal components of elements in V h ~ are 
linear. Each function v in the space V h ~ can be determined from its values at each interior 
vertex of the mesh, and of the mean value of the tangential component along each interior 
edge. In fact, in analogy with the discussion of the Mini element above, the space V h ~ can be 
decomposed as V h ~ = V h x © Vfr. As above the space V h l is the space of continuous piecewise 
linear vector fields, while the space in this case is spanned by quadratic "edge bubbles." 
To define this space of bubbles we let 

A 1 (T h )=A\(T h )uA 9 1 (T h ) 

be the set of the edges of the mesh Th, where A^^) are the interior edges and Af (Th) are 
the edges on the boundary of Q. Furthermore, if T G Th then A X (T) are the set of edges of 
T, and A\ (T) = A 1 (T) n A\ (T). 

For each e G A 1 (7/ l ) we let Q e be the associated macroelement consisting of the union of 
all T G Th with e G A±(T). The scalar function b e is the unique continuous and piecewise 
quadratic function on Q e which vanish on the boundary of fl e , and with J e b e ds = \e\, where 
|e| denotes the length of e. The space V% is defined as 

Vh = span{6 e t e | e G A\(T h ) }, 

where t e is a tangent vector along e with length |e|. Alternatively, if Xi and Xj are the vertices 
corresponding to the endpoints of e then the vector field ip e = b e t e is determined up to a 
sign as ip e = 6\i\j(xj — Xi), where {Xi} are the piecewise linear functions corresponding to 
the barycentric coordinates, i.e., Xi(xt) = <$i,fc for all vertices xj~- In particular, 

■0 e " ( x j — x^ ds = Q / XiXj ds\e\ 2 = \e\ 3 . 

J e 

As above the desired interpolation operator 17^ will be of the form 

n h = n b h (i-R h ) + R h , 

where Rh '■ L 2 (Q;M. 3 ) — > V h x is the same Clement operator as above, and where Il\ : 
L 2 (f2;IR n ) — > needs to be specified. In fact, to perform a construction similar to the 
one we did for the Mini element it will be sufficient to construct Il\ such that it is L 2 -stable, 



and satisfies the commuting relation (4.5). 



12 



KENT-ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER 



We will need to separate the triangles which have an edge on the boundary of from the 
interior triangles. With this purpose we define 

T h 9 = {TeT h \Tndne Af(T h )} and T h =T h \Tt 

In order to define the operator JJ\ we introduce Z h as the lowest order Nedelec space with 
respect to the mesh Th- Hence, if z G Zh then on any T G Th, z is a linear vector field such 
that z(x) ■ x is also linear. Furthermore, for each e G Ai(%) the tangential component of z 
is continuous. As a consequence, Z h C if (curl; Q) where the operator curl denotes the two 
dimensional analog of the curl-operator given by 

curlz = curl(*i, z 2 ) = d x ^z x - d Xl z 2 . 

It is well known that the proper degrees of freedom for the space Zh is the mean value of the 
tangential components of v, v • t, with respect to each edge in Ai(7^). Furthermore, we let 



z° h = {zez h \ [ z-tds = o,T eT h 9 }. 

J&T 



Alternatively, the elements of Z® are those vector fields in Zh with the property that 
curler = if T G T h 9 , i.e., z is a constant vector field on T for T G T h d . It is a key 
observation that the mesh assumption given above, that any T G Th intersects dfl in at most 
one edge, implies that the spaces Vfi and Z® have the same dimension. Furthermore, we note 
that gradg G Z® for any q G Qh- 

For each e G Ai(7h) let <p e G Z h be the basis function corresponding to the Whitney form, 
i.e., e satisfies 



(0 e • t e ) ds — |e|, and ^ (<j) e ■ t e >) ds = 0, e ^ e, 



where, as above, t e is a tangent vector of length e. Hence, if e = (x i? then the vector field 
e can be expressed in barycentric coordinates as 

e = A, grad Aj — Aj grad Aj. 

Any zeZj can be written uniquely on the form 

Z = ^ °e^e + ^ Ce< ^ e ' 
eeA\(T h ) e&Af(T h ) 

where the coefficients a e corresponding to interior edges can be chosen arbitrarily, but where 
the coefficients c e for each boundary edge should be chosen such that curl z = on the 
associated triangle in T h 9 ■ We note that there is a natural mapping §h between the spaces 
Vfi and Z® given by $/ l (^' e ) = 4> e for all interior edges, or alternatively, 



(4.8) J$ h {v)-t e ds = \e\- 2 lvt e ds, eeA[(T h ) 



Below, we will use Vfi(T) and Z^(T) to denote the restriction of the spaces Vfi and Z° to a 
single element T, and $ T will denote the corresponding restriction of the map & h . 

We will define the operator n\ : L 2 (Vt; W 1 ) ->• V% by, 
(4.9) (n b h u, z) = (u, z) \/zez° h . 
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To show that this operator is well-defined the following general formula for integration of 
products of barycentric coordinates over a triangle T will be useful (cf. for example [TJJ 
Section 2.13]) 

(4.10) / \?\?\?dx= m m 

where a\ = ai!a; 2 !a 3 !, \a\ = Y2i a i anc ^ 1^1 ^ s the area °f ^- 

Lemma 4.1. Lei T E T£ with edges ei,e 2 ,e 3 . For any v E V£(T) we have 

[ v ■ § T (v)dx > -laHll 
Jt 5 

where $ T = ^h\t, v = ^a*^ and \a\ 2 = X]i a i- 
Proof. A direct computation gives 

n 3 3 n 3 3/i 

/ v • $ T {v) dx = OjOj / ^ • ei ^ = «i% / b et (j) ej -t ei dx = a T Ma, 

Jt i=l j=l Jt i=l 3=1 Jt 

where the 3x3 matrix M is given by 

M = {M t Jl J=1 = {[ ^ • e , dar}? J=1 { / ■ t £i dx}{ J=1 , 

JT JT 

The desired result will follow from the diagonal dominance of this matrix. 

Let Xj be the vertex opposite ej, and let Xj be the corresponding barycentric coordinate 
on T. Then the first diagonal element M^i of the matrix M is given by 

M u = 6 / A 2 A 3 (A 2 gradA3 - A 3 grad A 2 )(a; 3 -x 2 )dx = 6 / (A 2 A 3 + A 2 A 3 ) dx = 2|T|/5, 
Jt Jt 



where we have used formula (4.10) in the final step. Actually, from this formula we derive 
that all the diagonal elements are given by M^i = 2|T|/5, and similar calculations for the 
off-diagonal elements gives |Mjj| = |T|/10. In addition, the matrix M is symmetric. The 
matrix M is therefore strictly diagonally dominant, and by the Gershgorin circle theorem 
all eigenvalues are bounded below by \T\/5. By combining this with the fact that M is 
symmetric we conclude that a 1 Ma > |a| 2 |T|/5, and this is the desired bound. □ 

The next lemma is a variant of the result above for T E T®. 

Lemma 4.2. Let T E T® with interior edges e±, e 2; and where e 3 is the edge on the boundary. 
For any v E V^{T) we have 

v ■ ® T (v)dx > ^\a\ 2 \T\, 
where v = aii[) ei + a 2 ip e2 and \a\ 2 = a\ + a\. 



14 



KENT-ANDRE MARDAL, JOACHIM SCHOBERL, AND RAGNAR WINTHER 



Proof. Let v = aiipe 1 +a2'4 , e2 — 6aiA2As(x3— X2)+6a2AiAs(x3— xi). It is a key observation that 
in this case $t(v) is simply given as $(u) = — aigradA2 — 02 grad Ai. Since grad A; • t ei = 
we therefore obtain from (4.10) that 

v ■ $ r (t>) dx = 6a 2 I A2A3 grad A2 • (#2 — £3) dx + 6a 2 , / A1A3 grad Ai • {x\ — x 3 ) dx 
Jt Jt 

= 6a 2 I A2A3 ax + 6a 2 , / A1A3 dx = -\a\ 2 \T\. 
Jt Jt 2 

This completes the proof. □ 
Lemma 4.3. There is a positive constant cq, independent of h, such that for each v 6 V£ 

(v,z) 



sup 

zez° h \\ z \\L2(n) 



> co\\v\\ L 2(n)- 



Proof. Let v G Vfi be given, i.e., v = Y^ e eA{(T h ) a ^ e - We srm ply choose the corresponding 
z = $h(v) = J2eeA\(T h ) a ^ + S eG Af (T h ) c ^ G Z h- lt follows from scaling and shape 
regularity that the two norms of z, given by 



(4.11) and (J^ Yl a ^ 



2^1/2 

TeT he eA[(T) 

are equivalent uniformly in h. Correspondingly, the two norms 

Te% eeA[(T) 



(4.12) IHmn) and (^ |T| 2 ^ 



are uniformly equivalent. As a consequence of these properties, combined with Lemmas |4.1 
and 4.2 , we obtain 

(v,$ h (v))=^2 I ® T {v) -vdx > - \ T \ Yl °^ - c o\\®h(v)\\ L 2 {n) \\v\\ L 2 (n) . 



Ten' 1 TeT h eeA\(T) 

where Co > is independent of h. This completes the proof. □ 



It is a direct consequence of Lemma 4.3 that the operators Il\ are uniformly bounded in 



£(L 2 (f2; IR 2 )). In fact, the associated operator norm is bounded by Cq . Note that in contrast 
to the situation for the Mini element, the operator 1T\ is not local in this case. However, if 
the mesh is quasi-uniform we obtain from (4.7) that 



(4.13) \\n b h u\\ HH n) <c{J2 h T 2 \\n^u\\h iT) ) 1/2 < crf^h-^nluW^a). 

Ten 



As a further consequence we now obtain. 



Theorem 4.4. The operator TI^ satisfies properties (4.3) and (4.4). 
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Proof. To show that the operator Ilh fulfills the condition (4.4), it is enough to show that 
the operator H\ satisfies the corresponding condition (4.5). However, since gradQ^ C Z\ 
this follows exactly as before, since 

(div II^v,q) = — (n^v , gradg) = — (t>,gradg) = (divt>,g). 



Furthermore, to show (4.3) it is enough to show that Ilh is uniformly bounded with respect 
to h in both £(L 2 (fi;lPj) and C(Hq(Q; IR 2 )). However, the L 2 -result follows from the 
corresponding bounds for the operators H h h and Rh- Finally, by combining (4.6), (4.13) and 
the boundedness of H\ in I? we obtain 

\\n h v\\ HHa ) < \\n b h {I - R h )v\\ HHn) + \\R h v\\ H i {Q) 

< cih^WnKl - R h )v\\ L 2 {n) + ||u|| H i(n)) 

< c(c 1 /i _1 ||(J - R h )v\\ L2{n) + \\v\\ H i in) ) 

< c||u||Hi(n), 

and this is the desired uniform bound in £(Hq(Q\ IR 2 )). □ 



4.1.3. More general triangulations. The analysis of the Taylor-Hood method given above 
leans havily on the assumption that there are no triangles in Th with more than one edge 
on the boundary of Q. This assumption simplifies the analysis, but it is not necessary. The 
purpose of this section is to relax this assumption. 

We let T^' 1 and T h 9 ' 2 denote the subset of triangles in Th with one or two edges in d£l, 
respectively. We let T h 9 = T h 9 ' 1 U T 9 ' 2 be the set of all boundary triangles, and as before 
T h l = T h \ T h 9 . We note that T G T h ' 2 then there is a unique associated triangle T G Th 
such that Tnr e A\(Th)- We will denote this interior edge associated any T E T h 9,2 
by ex, and we will use T* to denote the macroelement defined by the two triangles T and 
T~. The set of all interior edges of the form e^, T G T 9 ' 2 , will be denoted A z { 2 (Th), while 
A l { (Th) = A\{Th) \ A\ (%,)■ Throughout this section we will assume that all the triangles 
of the form T~ are interior triangles, i.e., 

T~ G T h l for all T G T 9 ' 2 . 

The interpolation operators Ilh and U\ will be defined as above, with the only exception 
that the definition of the bubble space V^ 6 is chang ed s lightly in the neighborhood of the 

i-2, 



edges in A^' (Th)- We note that the result of Lemma 



4.2 



the result of Lemma 4.3 a nd a s a consequence Theorem 



4.2 



still holds if T G T 9 ' 1 - To establish 
4.4, in the present case we basically 



for triangles T G T h ' • More precisely, we need to modify the 



need an anolog of Lemma 

definition of the map used in proof of Lemma [473] to the present case. Actually, the map 



will not be defined locally on the triangles T G T 9 ' 2 , but rather on the corresponding 
macroelement T*. As in the previous section the space Z® is taken to be the subspace of the 
Nedelec space Zh such that z G Z\ is constant on the triangles in T h 9 ■ To be able to define the 



interpolation operator mapping into the bubble space, by (4.9), the two spaces and 
Z® must be balanced. In particular, they should have the same dimension. However, in the 
present case the dimension of the space Z° is not the same as the number of interior edges. 
To see this just consider the restriction Z^(T*) of Z® to a macroelement macroelement T*, 
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Figure 2. The macroelement T* . 

cf. Figure [2] below. The dimension of the space Z®(T*) is four, while there are only three 
interior edges, namely the three edges of T - . To compensate for this we will extend the 
space of bubble functions V^, by including also "normal bubbles" on the edges e*r, T G 7^f' 2 . 

In the present case we define the space C Vh by 

V h h = span({ b e t e I e G A\(%) } U { b e n e | e G AfCft) }). 

Here t e and n e are tangent and normal vectors to the edge e with length |e|. With this 
definition the space V£ has the same dimension as Z%. Furthermore, the map : Vfi — > Z® 



will be defined to satisfy property (|4.8|) for all edges in (7h). Note that this specifies $^ 



on all triangles in 7^, except for the ones that belongs to the macroelements T* , T G 7^f' 2 . 
To complete the definition of $^ we need to specify its restriction to each macroelement T*. 

Consider a macroelement T* of the form given in Figure [2} Here the edge ct has endpoints 
denoted by x\ and X2, the third boundary vertex of T* is xo, while the single interior vertex 
of T* is £ 3 . We will use to denote the edge opposite x i} i = 1,2, of the triangle T, while 
e~ are the corresponding edges of the triangle T~ . We let $t : V^(T*) — > Z°(T*) be the 
restriction of §h to T*. To be compatible with the definition of outside the macroelements 



T* the map $t has to satisfy condition (4.8) on the edges e { , i.e., 
(4.14) / $ T (v)-t e -ds = \e;\- 2 / v-t er ds, i = 1,2, 



where £ e7 is a vector tangential to e i . As a basis for the space Z®(T*) we will use the 
functions 

07 = AigradA 3 - A 3 gradAi, i = 1, 2, 
with support only in T~ , combined with the two functions 0j given by 

grad A, on T, 
(-l)tyr, onT" 
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for i = 1,2, where 0t = Aj grad A2 — A2 grad Ai corresponds to the Whitney form associated 
the edge er = (xi, x 2 ). The functions 0j, 0~ for i — 1, 2 spans the space Z®(T*). 



We will define two basis functions ipi of Vjf as a multiple of the scalar bubble function b 
namely, 

(4.15) ipi = b eT Wi = 6AiA 2 Wj i = 1,2, 

where the vectors Wi will be chosen below. Furthermore, the functions ip~ , are given as 

ip~ = 6\i\ 3 (x 3 - Xi) + (3(-iy\ 1 \ 2 (x 2 - xi), i = 1, 2, 
where (3 = 6|T"|/(5|T| + 4|T~|). We note that < < 3/2. 

The functions for i = 1,2 span the space V%(T*), and we define = <pi 

and ^rlV'j - ) = 07- A. map of this form will satisfy the compatibility condition (4.14) by 
construction. The motivation for the choice of the constant is that we obtain 



(4.16) 



ipi ■ (j)jdx = 0, i,j = 1,2. 



Lemma 4.5. The orthogonality conditions (4.16) hold 



Proof. The identities (4.16) can be verified by formula (4.10). For example 



tp 1 ■ 0i dx = — (3 / AiA 2 grad Ai • (x 2 — Xi) dx 



T 



- J (Ai grad A 2 - A 2 grad Ai) • (6AiA 3 (a;3 - x x ) - /3AiA 2 (x 2 - x ± )) dx 



f3 / \t\2dx- / (6A1A2A3 - l3(Xf\ 2 + A1A2)) dx 
Jt Jt- 

(-6|r~| + /3(5|T| + 4|T _ |)/60 = 0. 



Furthermore, it is easy to check that 

ipi ■ 02 dx 



tp 1 ■ 0i dx, 

T* JT* 

and as a consequence § T * ipi ■ <p 2 dx = 0. Similar computations can be done for the integrals 
involving . □ 



We can also verify, again using formula (4.10), that the 2x2 matrix M = {M { ,}ij=i,2 



M~ 



IT" 



(24 - /3)/60 (6 + /3)/60 
(6 + /3)/60 (24 - /3)/60 



For < (3 < 3/2 this symmetric matrix is strictly diagonally dominant with both eigenvalues 
greater than |T _ |/4. 
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%3 



Figure 3. The macroelement T* and the associated parallelogram. 

Finally, we need to investigate the 2x2 matrix M = {Mij}ij = \p = {f T *ipi ■ 4>j}i,j=i,2- 
However, first we need to define the functions ipi = $t (</>«) precisely by specifying the vectors 



Wi in (4.15). We let 



■01 = 67AiA 2 (x 3 - x 2 ), and -02 = 67 1 A 1 A 2 (x 3 - x{), 
where the positive constant 7 will be chosen below. 

Assume for a moment that T* is a parallelogram. Then X 3 — X\ — X2 — Xq and x 3 - x 2 = 
x\ — xo, and therefore we would have easy computable representations of the functions ipi 
on both T and T~ . In general, we introduce a new point xq € IR 2 , depending on T", with 
the property that xq, x%, x 2 , x 3 corresponds to the corners of a parallelogram, cf. Figure [3j 
More precisely, 

£0 = x± — (x 3 - x 2 ) = x\ + x 2 - x 3 . 

Let {Aj} 2 =0 be the barycentric coordinates with respect to the triangle T, extended to linear 
functions on all of 1R 2 . Then Ai(x 3 ) + A 2 (x 3 ) > 1 and 



Xi{x ) + X 2 (x ) 



Ai(x 3 ) + A 2 (a; 3 ) < I. 



In fact, it is a consequence of shape regularity that there is a constant a > 0, independent 



of h and the choice of T e Tjf' 2 , such that 

(4.17) Ai(x ) + A 2 (x )<l-a. 

If we compute the matrix {f T ipi ■ 4>j dx} we obtain 



7(1 - Ai(x )) -7A2(£o) 
-T l k{x ) T\l-\ 2 (xv)) / 
To control the full matrix M we also need to consider the contributions from the triangle T~ . 



{J Ipi- (f)j dx}ij =lt2 



\T\ 
~2~ 



A straightforward computation, using formula (4.10), shows that the matrix {J T _ ipi ■ <pj} is 
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given by 

\T~\ I 7 



{ / ipi- <pj dx} iJ=h 2 



We will utilize the constant 7 to obtain a symmetric matrix M. We define 



7 



'2\T-\ + 5|T|Ai(£ ) 



2\T~\ + 5|T|A 2 (i ) 

This choice of 7 is motivated by the desired identity 

|T| - \T~ I |T| - IT" 

7 (^A 2 (x ) + = 7 _1 ( V A i( £ o) + 



2 v u/ 5 y ' v 2 v ' 5 n 
which can be seen to hold, and therefore the matrix M is symmetric. Furthermore, we note 
that 



7,7 1 < 




5]T| 
2\T- 



Therefore, it is a consequence of shape regularity that the positive constant 7 is bounded 

r 9 

'h 



1 2 

from above and below, independently of h and the choice of T e 7^ ' 



Lemma 4.6. T/ie matrix M defined above is symmetric and positive definite with both eigen- 
values bounded below by C\\T\, where C\ = amin(7, 7 _1 )/2. 



Proof. It follows from the calculations above that 
M = 



7(^(1 - X 1 (x )) + l» -7(^A 2 (x ) + ^ 



^(^^(xo) + V) 7^(^(1 - A 2 (x )) + » 

Since Ai(xo) + A 2 (xo) < 1 — a it follows from Gershgorin circle theorem that both eigenvalues 
of M are bounded below by ^p- min(7, 7 _1 ). □ 

We now have the following result. 
Lemma 4.7. The conclusion of Lemma \4 ■ 3*| holds in the present case. 

Proof. Let v G be given. We first consider the situation on each macroelement T*. If 
v = 'Yjii. a t' l Pi + a ii ) ^) e Vh(T*) then we write v — v~ + v + , where vT = J7 a^ip~ . Observe 



that Lemma 4.6, together with the orthogonality property (4.16), implies that 

/ v-$ T (v + )dx= v + ■ $ T {v + )dx > c x \T\\a + \ 2 . 
Jt* Jt* 

Similarly, we have from the property of the matrix M~ , the norm equivalences expressed by 



(4.11) and (4.12), and shape regularity that 

/ v ■ $ T (v-) dx > / v- ■$ T (v~)dx-\\v + \\ L 2 {T - ) \\$ T (v-)\\ L 2 (T - ) 
Jt- Jt- 
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> t IT - 1 la" 

> \\T~\\a 



c\T \\a \\a 

IT* I I +|2 

c 2 \T a , 



where the constant C2 is independent of h and T. By choosing <&t{v) = C$(t> + ) + $(u 
where the constant C is sufficiently large, we can now conclude that 

v ■ $ T (v) dx > c\T*\\a\ 2 . 



We note that the map <&t will inherit the compatibility condition (4.14) from the map $t 



By combining this result on each macroelement T*, with the map $^ defined previously on 
the rest of the triangles in Th, to a global map <i\ mapping Vfi, we can conclude, as in the 



proof of Lemma 4.3 that 



(v,z) 

SUP j—r 

zez° \\ z \\L2(n) 



> 



(v,$h(v)) 
\$h(v) 



> C \\v\\ L 2 (n) . 



This completes the proof. 



□ 



As we have noted above the result just given implies that the conclusion of Theorem AA_ 
holds will hold for the more general meshes studied in this section. 
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